Topological order in a three-dimensional toric code at finite temperature 
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We study topological order in a toric code in three spatial dimensions, or a 3+ ID Z2 gauge theory, at finite 
temperature. We compute exactly the topological entropy of the system, and show that it drops, for any in- 
finitesimal temperature, to half its value at zero temperature. The remaining half of the entropy stays constant 
up to a critical temperature Tc, dropping to zero above Tc. These results show that topologically ordered phases 
exist at finite temperatures, and we give a simple interpretation of the order in terms of fluctuating strings and 
membranes, and how thermally induced point defects affect these extended structures. Finally, we discuss the 
nature of the topological order at finite temperature, and its quantum and classical aspects. 



I. INTRODUCTION 

Some quantum systems are characterized by a type of order 
which cannot be captured by a local order parameter that sig- 
nals broken symmetries, but instead the order is topological 
in nature. 1 1] One of the ways in which this topological order 
manifests itself is in a ground state degeneracy that cannot be 
lifted by any local perturbation, and that depends on the genus 
of the surface in which the system is defined. Recently, there 
have been efforts to find characterizations of topological order 
other than ground state degeneracies, in particular exploring 
the entanglement in the ground state wavefunction. 

At zero temperature, topological order can be detected us- 
ing the von Neumann entanglement entropy, more precisely 
a topological contribution to it that can be separated from the 
boundary contribution by appropriate subtractions of different 
bipartitions of the system. Because the pure state density 
matrix is constructed from the ground state, it was argued in 
Ref.|2 that topological order is a property of the wavefunction, 
and not of the Hamiltonian, at absolute zero temperature. 

An interesting question is what happens with topological 
order at finite temperature. The question is relevant because 
thermal fluctuations, no matter how small, are present in any 
laboratory system. To address this issue, it was proposed in 
Ref. 1 4] to use the topological entropy as a probe of topolog- 
ical order, but to compute it using an equilibrium mixed state 
density matrix p = Z^'e^P^. It becomes clear that, as op- 
posed to zero temperature for which one can do away with 
the full information contained in the Hamiltonian and just use 
the ground state wavefunction, topological order, if present at 
finite temperature, must be a property of the Hamiltonian. 

The topological entropy was computed exactly for the 2D 
Kitaev model at finite temperature T, and it was shown 
that the infinite system size limit and the T ^ limit do not 
commute, and that at finite T the topological entropy vanishes 
in the thermodynamic limit. Thus, it was argued that the topo- 
logical order in the 2D system was fragile. [4i,|7|,[8(] 

Here we show that the situation in 3D is rather different, 
using the 3D version of Kitaev's model as an example. |9] 
In contrast to 2D, topological order survives up to a phase 
transition at a finite temperature T^. The order can be probed 
through a non-vanishing topological entropy, as well as un- 
derstood from a simple cartoon picture that we present in the 



paper, using the fact that in 3D strings can move around point 
defects (as opposed to 2D). 

We prove in this paper that the von Neumann entropy of a 
subsystem j4. of a Z2 gauge model such as Kitaev's toric code, 
in any number of dimensions, can be always decomposed into 
two additive contributions from each of the two gauge struc- 
tures (magnetic and electric): idoll 

SyMT) = S^^'l{A;T/XA) + si,''^{A;T/Xs), (1.1) 

(s) (P) 

where ^yj^, and ^yj,, are the separable contributions from the 
stars and plaquettes of the model, and Xa and Xb the associated 
coupling constants for these two structures. Consequently, the 
same additive separability holds for the topological entropy, 
which is a sum of two independent contributions: 
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One of the contributions, S^^p^, evaporates for any infinitesi- 
mal temperature in the thermodynamic limit, just as in 2D, but 

(P) 

the other one, Sf^^^, remains constant up to a finite tempera- 
ture phase transition at = 1.313346(3)A-b, that occurs for 
the 3D case: 



n3D 



21n2 T = 

,„p„(r) = <! in2 < r < 

T> Tc. 



(1.3) 



As a consequence of these results, we argue that topologi- 
cal order can be well defined at finite temperatures in 3D. ifllll 
This finding raises the following interesting question: is the 
finite T order classical or quantum? Perhaps another way to 
ask the question is the following: Which kind of information 
can be robustly stored using the isolated topological sectors 
in phase space that cannot be connected by local moves (2^ 
such states in 3D): classical (bits) or quantum (qubits) infor- 
mation? While we cannot argue that the system does not real- 
ize a full quantum memory, we can at the least argue that it can 
store probabilistic information (pbits - probabilistic bits ifTsll ) 
in the form of a quantum superposition of states in the dif- 
ferent topological sectors, where the square amplitudes for all 
states in a given sector (a probability) does not fluctuate in the 
thermodynamic limit if the coupling to a thermal bath is local. 
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However, the relative phases for all these amplitudes could be 
scrambled. This weak type of quantum superposition is not 
discernible from a classical probability distribution. 

Finally, this example shows that the notion of classical 
topological order, suggested for hard constrained models in 
2D, IMI is well defined in 3D without resorting to any hard 
constraints. 



II. THE MODEL 

Consider a three-dimensional version of Kitaev's toric 
code, ^ defined on a simple cubic lattice of size N ~ Lx Lx 
L, with periodic boundary conditions and spin- 1/2 degrees of 
freedom a, living on the bonds, ; = 1 , . . . , 3A^ (a^, oj and of 
being the three Pauli matrices). Let us label the centers of 
each single square plaquette in the lattice with p = 1 , . . . , 3A^, 
and each site of the cubic lattice with s = 1,...,N. 

Let us define the plaquette and star operators on the lattice 



(2.1) 



as illustrated in Fig.[T] The Hamiltonian of the model can then 




FIG. 1: (Color online) - Illustration of the Kitaev model in 3D, with 
explicit examples of a star operator A, = O*^/ at the lattice site s, 
and of three plaquette operators Bp — ]Ja^ at the plaquette-dual sites 
Pi , P2 and p3. The a spin index ; labels respectively the 6 (red) spins 
around s and the 4 (blue) spins around p (connected by dashed lines). 



cube, so we have one less constraint). Moreover, three addi- 
tional constraints come from the fact that the product of all 
plaquette operators along any crystal plane in the cubic lat- 
tice (i.e., {x,y), {x,z), or {y,z)) yields the identity, and we are 
finally left with 2N — 2 independent plaquette operators. 

The ground state (GS) manifold of the system is identified 
by having all plaquette and star quantum numbers equal to 
+ 1, and it is 2^'^^''^^')^^^^^) — 2^ dimensional, assuming 
periodic boundary conditions in all three directions. Similarly 
to the 2D case, one can notice that this degeneracy has a topo- 
logical nature, and the different sectors are distinguished by 
three non-local operators 



or 



^1 



(2.3) 



n< (2.4) 



that are diagonal in the o' and cr^ basis, respectively. Here 
the Yj- can be any winding paths along the edges of the cubic 
lattice in each of the three crystal directions (x, y, or z), and 
the t,j can be any winding planes perpendicular to each of the 
crystal directions and passing through the midpoints of the 
corresponding edges of the cubic lattice (i.e., crystal planes in 
the dual lattice whose sites sit at the centers of the elementary 
cubic cells). Two examples are shown in Fig.|2]for clarity. 




be written in terms of these operators as 



(2.2) 



FIG. 2: (Color online) - Two examples of the non-local operators 
needed to distinguish between the degenerate GS of the 3D Kitaev 
model. 



where A,^ and are two real, positive constants. 

Notice that all star and plaquette operators commute, but 
they are not all independent. While only the product of all 
star operators equals the identity, therefore leaving — 1 in- 
dependent star operators, the product of the plaquette opera- 
tors around each cubic unit cell gives the identity, therefore 
introducing — 1 constraints in the 3A^ total plaquette opera- 
tors (the product of all but one cube is equivalent to that same 



In the & basis and in the topological sector where all the 
equal + 1 , the GS wavefunction of the system can be written 

as 



\GS) 



1 



10), 



(2.5) 



where |0) is any state in the sector, say the state with all the 
a~- ~ +1, and G is the Abelian group generated by all products 
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of star operators (of dimension \G\ = 2^~^). In the a* basis 
and in the topological sector where all the Z- equal +1, the 
GS wavefunction of the system can be written as in Eq. ( 12.5b . 
where now |0) is any state in the sector, say the state with 
all the a] — +1, and G is the Abelian group generated by all 
products of plaquette operators (of dimension \G\= 2^'^^''). 

Notice the two different underlying structures in the system: 
the closed a"' loops along the edges of the cubic lattice, which 
satisfy Hioop ^ ^ identically, and the closed O'^ membranes 
in the body-centered dual lattice (locally perpendicular to the 
edges of the original lattice), satisfying Hmembrane — ^ iden- 
tically (see Fig.O. 




versa - e.g., a spherical shell. Thus, there is no unique way 
to generalize the 2D case. Two equally valid options are il- 
lustrated in Fig.m based on a 'spherical' (1-4) and a 'donut- 
shaped' (5-8) bipartition scheme, respectively. 

In the O' basis |(5J, where G is generated by the star opera- 
tors, the calculation of the entanglement entropy Sy^ proceeds 
as in the 2D case, 10, [l3> [Htl • Using the group property of G 
in Eq. i2.5i . one can show that 



i-A) — 



(3.1) 



where dj^ is the dimension of the subgroup Gj^ C G con- 
taining all the elements of G that act as the identity on "B, 
Gj^ = {g ^ G \ g = gji^ Irg}, and similarly for subsystem S. 
As in the 2D case, these subgroup dimensions depend on the 

(s) (s) 

number ) of star operators acting solely on spins in 

A CB), and on the number m^j (m^) of connected components 
oiACB): 



,-1 



d^ = 2^^'B 



(3.2) 
(3.3) 



The contribution to c/^^. 



and vice versa the m^^ contribu 



tion to drg, come from the so-called collective operations, i.e., 
elements of the groups (G^) that cannot be expressed as 
products of star operators in A (25). In the 3D case, such col- 
lective operations correspond to non-contractible closed mem- 
branes. In this respect, bipartitions 1 and 8 are special in that 
subsystems 123 and 8^1 are composed of two separate con- 
nected components {m^^ — m^j^ = 2), while all other subsys- 
tems have only one component. 

We can then compute the topological entropy 5(opo of the 
system in the basis from either the spherical or the donut- 
shaped bipartition scheme. 
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FIG. 3: (Color online) - Two examples of the underlying struc- 
tures of the 3D Kitaev model: the closed CT~ loops along the edges of 
the cubic lattice, which satisfy Oioop^; ~ ^' closed cr*^ mem- 

branes in the body-centered dual lattice, satisfying rTmembrane ^/ ~ ^■ 



where we used the fact that all IN'^*' contributions cancel out 

exactly. In fact, if we define J^f^s = - '^il - ^° be 
the number of star operators acting simultaneously on A and 
■B, J^^"^ = N being the total number of star operators in the 
system, one can show that 



III. THE TOPOLOGICAL ENTROPY AT ZERO 
TEMPERATURE 

Let us first compute the zero-temperature topological en- 
tropy of the system, using a three-dimensional version of the 
bipartition scheme proposed by Levin and Wen in two di- 
mensions. Notice, however, that in 3D a bipartition can be 
topologically non-trivial with respect to closed loops but not 
with respect to closed membranes - e.g., a donut -, and vice 
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(3.5) 



This result relies on the fact that the total boundary in bi- 
partitions 1 and 4 is the same - with the same multiplic- 
ity, and with precisely the same edge and corner structure 
- as in bipartitions 2 and 3, by construction. Therefore, 



(1) 



(2) 



(3) 



(4) 




FIG. 4: Illustration of the two bipartition schemes used for the 3D Kitaev model: 'spherical' (Top), and 'donut-shaped' (Bottom). 



Similarly for bipartitions 



5-8. 

Let us also compute the topological entropy in the ba- 
sis, as it will be useful when we consider the finite tem- 
perature case. The group G is now generated by the plaquette 
operators, which are highly redundant and require more in- 
volved calculations to obtain the von Neumann entropy Syj^. 
In fact, while Eq. (13. Il l still holds, one needs to count the num- 
ber of independent plaquette generators of subgroups Gj^ and 
G3 in order to obtain the equivalent of Eqs. (I3.2l i and ( I3.3l l. 
Notice that the collective operations are now given by closed 
loops, and only bipartitions 4 and 5 allow for non-trivial (i.e., 
non-contractible) loops. 

As we discussed before, \G\ — 2^'^^''. This arises from 
counting all independent generators of G as the total number 
of plaquettes in G (all possible generators), minus the number 
of independent constraints. These are all but one of the cubic 
unit cells, plus three crystal planes. Similar arguments apply 
to the bipartitions 1-8. Notice that in all of the bipartitions, 
subsystem A does not contain any entire crystal plane, while 
subsystem S always contains all three crystal planes. Taking 
advantage of this simplification, in the following it will be 
understood that G3 has three less independent generators with 
respect to . 

Let us proceed case by case. For bipartitions where both 
A and S have only one connected component without han- 
dles, such as bipartitions 2,3,6, and 7 in Fig.|4] the group Gy^ 
(equivalently Grg) is generated by all the plaquette operators 
acting solely on A, subject to the constraints given by all cu- 
bic unit cells entirely contained in A. There are no collective 
operations in this case, and one obtains 



A 



2^33 



(3.6) 
(3.7) 



where DnT^^ is the number of plaquette operators acting on 



spins in A, JvT^^ is the number of cubic unit cells in A, and 
similarly for S. 



Consider then the case of bipartition 4 (equivalently, 5). 
Although both A and 23 are still connected, the presence of 
a handle allows now for collective operations. Take a crys- 
tal plane perpendicular to the largest surface of subsystem A, 
and draw it so that it bisects the donut into two identical U- 
shaped portions [see Fig. |5] (Top)]. The intersection of this 




FIG. 5: (Color online) - Illustration of the collective operations in 
bipartitions 4 and 5 in the a' basis, acting on subsystem 3 (Top) and 
subsystem A (Bottom), respectively. 

plane with A gives two rectangles of size rx{R — 2r), a dis- 
tance R — 2r apart. Now take the product of all plaquettes 
belonging to one of the rectangles plus those at its boundary. 
The resulting operation acts on 15 alone, yet it cannot be con- 
structed from plaquettes in S because the "outer boundary" 
of the rectangle cannot be the sole boundary of a surface in 
23. Notice that this collective operation can be deformed at 
will and moved along the donut by appropriate products of 
plaquettes in 23, therefore there is only one independent such 
operation. Similar arguments apply if we repeat the construe- 
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tion starting from a plane parallel to the largest surface of the 
subsystem A, again chosen so as to bisect the donut. This 
yields another independent collective operation acting now on 
A [see Fig.|5](Bottom)]. As a result, 



(3.8) 
(3.9) 



where nj^ = \ and n,^ = 1 are the number of collective opera- 
tions in A and S, respectively. 

Finally, one can show that there are no collective opera- 
tions in the cr*^ basis in bipartitions 1 and 8. In fact, all closed 
loops are contractible to a point both in A and in "B in these 
bipartitions. However, the disconnected nature of subsystem 
23 in bipartition 1 (equivalently, subsystem A in bipartition 8), 
requires special care in the counting of the independent gen- 
erators of Gj^ (respectively, G^). As in the previous cases, 
all plaquettes in A belong to Gy^, and all cubic unit cells in A 
act as independent constraints towards the counting of the in- 
dependent generators of Gj^ . However, in bipartition 1 , there 
is a class of closed membranes in A that cannot be assem- 
bled as a product of cubic cells in A. This is the case, for 
example, of the closed cubic membranes in A that surround 
entirely the inner component of 25. Any two such membranes 
can be obtained one from the other via multiplication by cu- 
bic unit cells in A. Thus, they only give rise to one additional 
constraint in the counting of the independent generators. In 
general, the number of such constraints is given by — 1, 
where is the number of connected components of 23. Sim- 
ilarly for bipartition 8 and subsystem 23, one obtains my^ — 1 
additional constraints, where nij^ is the number of connected 
components of A. 



Combining all of the above considerations into a general 
expression for the dimensions of subgroups Gy[ and in 
the CT^ basis, one obtains 



2^B 



(p) 



(3.10) 
(3.11) 



where ?»^''''' (m^ '''') is the number of crystal planes (c.p.) 
entirely contained in A (23). Recall that all bipartitions of in- 
terest have wi^'^ ' = and m^ '' ' — 3. 

We can then use Eq. ( 13.11 ) to compute the topological en- 
tropy of the system using the spherical and the donut-shaped 
bipartition schemes in the a^' basis. 
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(3.12) 



where we used the fact that all "N^p^ and M^'^' contributions 
cancel out exactly. In fact, if we define = J^^''' 



J^l^ to be the number of plaquette operators acting simul- 
taneously on A and 23, Tsf^''' — 3N being the total number 
of plaquette operators in the system, and we define ^Nfj^'^ = 
j^U') — j\f|^' _ jsj'l^'' to be the number of cubic unit cells simul- 
taneously encompassing spins in A and in 23, N^'' = being 
the total number of cubic unit cells in the system, one can 
show that 



jrip) _ jric) 



■>^4A ■'^4A 
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This result relies on the fact that the total boundary in biparti- 
tions 1 and 4 is the same - with the same multiplicity, and with 
precisely the same edge and corner structure - as in biparti- 
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(P) _ 



iAT,- 



Simi- 



tions 2 and 3, by construction. Therefore, "N 

^iL + and + = ^\ 

larly for bipartitions 5-8. 

Clearly, both bipartition schemes capture the topological 
nature of the system, and provide an equally valid measure 
of the topological entropy. In 2D the choice of bipartitions 1- 
4 in Ref. 2 is such that bipartition 1 is topologically equivalent 
to bipartition 4 upon exchange of subsystem A with subsys- 
tem 23, while bipartitions 2 and 3 are actually topologically 



invariant upon the same exchange. Hence, because the von 
Neumann entropy for the ground state is symmetric under the 
exchange of A and 23, the topological contribution measured 
in the 2D scheme is bound to be double counted, namely 
•Sfnnr, = 21nD = InD^, where D is the so called quantum di- 
mension of the system. |12|,|3I] In 3D, both the scheme 1-4 and 
the scheme 5-8 isolate the topological contribution to the en- 
tanglement entropy without double counting. Notice that all 
the bipartitions are topologically invariant under the exchange 
of A and 23, except for bipartitions 1 and 8. If we want to 
recover the symmetry of the 2D scheme, a possible solution is 
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to define 
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with D — 2. As we will see in the following, the symmetric 
1-8 choice is actually required if we are interested in studying 
the finite temperature case, since the von Neumann entropy 
is no longer invariant upon exchange of A and S, and a non- 
topologically-symmetric choice of bipartitions would lead to 
different results depending on whether we work with subsys- 
tem A or subsystem "B. ||17[1 



IV. THE FINITE TEMPERATURE BEHAVIOR 

In this section we study the behavior of the entanglement 
and topological entropies at finite temperature, via a general- 
ization of the approach used for the 2D Kitaev model in Ref. 4. 

A qualitative picture of the effect of thermal fluctuations 
can be argued by comparison with the two dimensional case. 
There the information about the topological sectors is stored 
in the eigenvalues of winding loop operators, namely products 
of spin operators along winding loops. On a torus, there are in- 
finitely many choices for such winding loop operators, but the 
absence of magnetic and electric charges (i.e., plaquettes and 
stars with eigenvalue — 1 ) in the gauge structure at zero tem- 
perature reduces them to only two independent ones: the two 
non-contractible winding loops on the torus. Any other can be 
obtained from these two via multiplication by an appropriate 
set of plaquette or star operators, which have eigenvalue +1 
at r = 0. Clearly the presence of order 1 (deconfined) ther- 
mal defects destroys immediately all topological information 
stored in the system, since the eigenvalues of two loops on 
opposite sides of a defect are no longer consistent with each 
other (see Fig.|6l). 

Let us now consider the case of the Kitaev model in 3D. 
First of aU, we need to discuss the two gauge structures sepa- 
rately, since they are no longer identical as in 2D. If we work 
in the a" basis, then the topological information is stored in 
the eigenvalues of winding membrane operators, given by the 
product of all operators belonging to a closed winding sur- 
face locally perpendicular to the bonds of the sites it crosses 
(see Fig. |2]i. All possible choices of these membranes yield 
the same result at zero temperature since the corresponding 
operators can be obtained one from the other by products of 
sets of star operators, which have all eigenvalue +1 in the GS. 
Thermal defects in this case play exactly the same role as in 
2D, since two membranes on opposite sides of a defect read 
off opposite eigenvalues of the corresponding winding mem- 
brane operator 

On the other hand, the situation is quite different for the 
loop operators defined in the O' basis. There the topological 
information is stored in winding loop operators - as in the 2D 
case - but they are now embedded in 3D. Clearly, localized de- 
fects have no disruptive effects on the topological information 




FIG. 6: (Color online) - Qualitative illustration of the disruptive ef- 
fect of two defects (solid red dots) in the 2D Kitaev model on a torus: 
two winding loops (black wavy lines) on either side of a defect (solid 
circle) read off opposite eigenvalues of the corresponding winding 
loop operator. 



because any two winding loops (with equal winding numbers) 
can be smoothly deformed one into the other without crossing 
any defects at low enough temperatures (see Fig. |2l). This is 




FIG. 7: (Color online) - Qualitative illustration of the reason why 
the topological information stored in the underlying loop struc- 
ture of the 3D Kitaev model is robust to thermal fluctuations: even in 
presence of sparse defects (solid red circles), any two winding loops 
(black wavy lines), with equal winding numbers, can be smoothly 
deformed one into the other without crossing any defects. (The wig- 
gly lines represent qualitatively the confining strings between defect 
'pairs' discussed in the text.) 

indeed the case here, where we learn from 3D lattice gauge 
theory that defective plaquettes are confined at low temper- 
atures. They are created in quadruplets by a single spin flip 
operation, and they can be pairwise separated only at the cost 
of creating a string of defective plaquettes in between the two 
pairs. 1 18, 19] Therefore, the winding loop operators will keep 
carrying the same quantum information in presence of a low 
density of defects. If we were to read out the topological infor- 
mation from the system, we would be getting the correct result 
as long as the chosen loop does not pass directly through a de- 
fect. 
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However, can this information be accessed by means of the 
same expectation values of loop operators that are used at zero 
temperature, Eq. (12.3b and Eq. ( I2.4l i? The answer to this ques- 
tion is negative, as it was recently shown using gauge the- 
ory arguments in Ref. J_2, A simple reason as to why naively 
choosing a given loop operator and looking at its expectation 
value alone does not capture the order below Tc is that, typi- 
cally, winding loops will pass through at least one defect in the 
thermodynamic limit (the probability of a loop not crossing 
any defect scales as (1 — Pdef)^^ where p^ef is the equilibrium 
density of defects at a given temperature, and L is the linear 
size of the system). However, only those loops that avoid the 
defects contain the topological information. (Recall that in 2D 
the eigenvalues of loop operators, even when they do not pass 
through defects, differ on two sides of one defect, in contrast 
to the situation in 3D.) This implies that the average expecta- 
tion value of loop operators is bound to vanish exponentially 
in system size for any finite density of defects, i.e., for any 
finite temperature, independently of the nature of the system. 
As we show in the following, the topological entropy of the 
system is capable of capturing these physical differences, and 
it accurately reflects the topological properties of the different 
phases. 

The physical meaning of the distinct sectors can be under- 
stood as follows. Consider preparing the system in a coherent 
superposition of different topological sectors at zero temper- 
ature. Raise the temperature to some value T < Tc, and then 
lower it again back to zero. If defects are confined, transitions 
between different loop sectors are forbidden throughout the 
process. We are thus bound to obtain a final state where the 
probability (magnitude of amplitude square) of finding the fi- 
nal state in each loop sector is the same as in the initial state. 
In this sense, the loop sectors are protected from thermal fluc- 
tuations at low temperatures, and topological order survives at 
finite temperature (T < Tc). 

That the system does not change sectors during the time that 
it is in thermal equilibrium with the bath is a dynamical prob- 
lem (broken ergodicity). This can be understood by contrast- 
ing the time scales for mixing sectors if defects are confined or 
deconfined. Deconfined thermal defects are free to randomly 
walk across the system, and induce transitions between differ- 
ent topological sector by means of creation, system-spanning 
propagation, and annihilation processes. The characteristic 
time for a sector-changing process scales therefore as some 
power of the system size, Xdeconfined ^ i"- In contrast, con- 
fined defects will have to overcome an energy barrier of the 
order of L to be able to wind around the system and induce 
a change in the topological sector. As a result, their char- 
acteristic time scale is instead exponential in system size, 
^confined ^ s'^- Even for rather small systems, confined defects 
would require time scales larger than the age of the universe 
to transition between sectors. 

An even more interesting situation occurs when both Z2 
gauge defect types are confined, so that the and topolog- 
ical sectors are both protected. This case is briefly discussed 
in Appendix lAl and it is related to error recovery that was ar- 
gued to be realizable for example in a 4D toric code. [Tj] What 
we argue here based on the finite temperature studies is that 



the system can be self-correcting: if the system is prepared in 
a given superposition at zero temperature and its temperature 
is raised and again lowered to zero without ever going above 
Tc, the system returns to the same original quantum state (a 
"boomerang" effect). 

The protection holds at low temperatures, but it is bound 
to vanquish as the density of defective plaquettes with eigen- 
value — 1 grows with temperature: once enough defects are in 
place, one can no longer deform paths around them. There- 
fore, we expect a loss of topological information as temper- 
ature is increased, via a topological phase transition at finite 
temperature. 

In analogy with 3D lattice gauge theory, we expect this 
transition to occur when plaquette defects deconfine at high 
enough temperature. This is captured by the expectation value 
of Wilson loop operators, which is exponentially suppressed 
with the length of the loop (perimeter law) at low tempera- 
tures, while it is suppressed with the area of the minimal en- 
closed surface (area law) at high temperatures. 1 18, 19, 20, 2ll 
In our notation, the transition temperature is set by the energy 
scale Xb, and the transition is expected to occur at the critical 
point of the 3D lattice gauge theory. 

The topological entropy is a non-local order parameter that 
detects the presence of topological order in a system. Any loss 
of topological information, e.g., whenever some topological 
sectors become ill-defined, should have a measurable effect 
on such entropy. Indeed, we show below that this is the case, 
and that the qualitative picture inferred from the arguments 
above is confirmed by an exact calculation of the topological 
entropy at finite temperature. 

A. The density matrix 

Let us work for convenience in the (tensor product) ba- 
sis, where the Hilbert space 3i is spanned by the whole set of 
orthonormal states |a), labeled by the configurations a of a 
classical Ising model on the bonds of a 3D simple cubic lat- 
tice (the value ±1 of each Ising variable corresponds to the 
eigenvalue of the operator at the same site). 

Define G to be the group generated by all plaquette opera- 
tors Bp = Yliep'^r Recall that any two elements of the group 
differing by products of plaquettes around closed membranes 
are in fact the same element (i.e., they are defined modulo the 
identities Flciosed membrane = 1)' where we are assuming pe- 
riodic boundary conditions, and full crystal planes are there- 
fore closed membranes as well. Recall also that \G\ = iP'^'^, 
where is the number of sites in the simple cubic lattice. Ev- 
ery two elements of the group commute with each other, and 
g^ = \, \/g £ G. For later convenience, let us label with a = 
the fully magnetized state o'^ = + 1 . 

The equilibrium properties of the system at finite tempera- 
ture are captured by the density matrix 

^ I„p(ak-P^|P) |a)(P| 
I„(ak-P^|a) • 
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For convenience of notation, let us rewrite the Hamilto- 
nian i2.2i as 

H = — — 

p = L^v 

p 

S - EA,. 

Notice that 5|a) = M^{a)\a), where M^(a) is the net "star 
magnetization", i.e., the difference between the number of 
stars with eigenvalue + 1 and with eigenvalue — 1 in the state 
I a). The action of any group element g is to flip plaquettes, 
which cannot change the sign of any star operator since they 
commute, and therefore Sg\a) — M^,(a)^|a), \/g e G. 

Thus, the denominator of Eq. ( 14. Il l becomes 

Y^{a\e-^"\a) = Y^e^'>^^''^'^^\a\e^'^BP\a). (4.2) 

a a 
Upon expanding 

gp^gf ^ j-j ^^Qgjj pj^^ _^ gjjjjj pj^^ ^ (4 3) 

p 

as follows from the definition P = Y.pBp and from the fact that 
B^-, = 1, one can explicitly compute the last term 

(a|eP^s^|a) = (ain [coshpA-B + sinhpA-gB,,] |a). (4.4) 

p 

All non-vanishing contributions in Eq. (14.4b come from 
products of plaquette operators that reduce to the identity (i.e., 
products around closed membranes). The above equation is 
therefore independent of a, which we set to the reference state 
in the following. 

The set of all possible closed membranes in a periodic 3D 
simple cubic lattice is in one-to-two correspondence with all 
possible configurations of an Ising model on the dual simple 
cubic lattice (the membranes are, say, the antiferromagnetic 
domain boundaries), provided we allow for both periodic and 
antiperiodic boundary conditions in all three directions. In 
this language, the sum of all non-vanishing contributions can 
be written as 

{Q\e^^BP\Q) = [coshpX5]-^'^-^AF(e) [sinhpA^l'^APle) ^ 

(4.5) 

where 6 is a generic configuration of the 3D Ising model with 
any type of boundary conditions, 3A^ is the total number of 
nearest-neighbor (nn) bonds and Np^{&) is the number of an- 
tiferromagnetic nn bonds. The factor 1/2 comes from the 
symmetry: a given membrane configuration corresponds to 
two equivalent but distinct Ising configurations. For conve- 
nience, let us introduce the simplified notation c = cosh pA-^, 
s = sinhpAg and t — s/c = tanhpA^, and define 7 > such 
that e^^'' = t (recall that > 0). The above expression can 



then be further simplified to 

e 




= {scr/'j:e^p(jj:sA 

e V {i-j) ) 
EE (scf^l^Zf, (4.6) 

where zy is the partition function of an Ising model on a sim- 
ple cubic lattice of size N — Lx Ly. L with reduced ferromag- 
netic coupling constant J, summed over all possible choices 
of (periodic or antiperiodic) boundary conditions. 1 16] 

We can now move on to compute the numerator of Eq. (14. 11 1. 

E(ak-P«|P)|a)(P| = 

= ][^eP^AM.(«)(a|eP^«^|p) I a) (PI 
a,p 

^Y^Y^e^W'^){a.\e^^BPg\a) |a)(a|^, (4.7) 

geG a 

where we used the fact that all matrix elements (a|eP^B^|p) 
vanish identically unless |p) = ^|a), 3^ G G. Once again, the 
expectation value {a\e^'^BP g\<x) is independent of a, and the 
above expression simplifies to 

££eP^AM.(«)(0|eP^B^g|0) \a){a\g. (4.8) 
geG a 

The expectation value can be computed explicitly by ex- 
panding the exponential 

= (0| n [cosh PA5 + sinh pA^ B J H !«) • 

P p'eg 

(4.9) 

Here, the notation Wpi^^B^, represents the decomposition of 
g in terms of the group generators {B^}. Clearly this de- 
composition is highly non-unique, since the group elements 
are defined modulo the identities noosed membrane 

Bp = 1, and 

Eq. ( 14. 9t needs to be handled with care. 

As before, all non vanishing contributions come from prod- 
ucts of plaquette operators that reduce to the identity. In this 
case, however, there are two options for every operator B^, : (i) 

it can be multiplied out directly by sinh PA^ B^, with p — p' 
(recall that B^y — 1); or (ii) it can be completed to an identity 
by an appropriate product of Bp terms so that B^, WB^ forms 
a closed membrane. Notice that in the second case the product 
over p may not include p' itself. 
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All this can be expressed in more elegant terms in the 
Ising language defined previously. Case (i) corresponds to 
the two spins across the bond being ferromagnetically 
aligned in the Ising model, and contributing a Boltzmann fac- 
tor sinh pA,^. Case (ii) corresponds to the two spins across p' 
being antiferromagnetically aligned, and contributing a Boltz- 
mann factor coshpA,^. Notice that the correlations between 
the different are automatically taken care of in the Ising 
language, and we obtain 

2 (OI.P^^^^IO) = (.c)^^/2^ exp (j £ n,j{g)StSj] 

EE (5c)3^/2zr(g), (4.10) 

where 

Recall that a bond in the Ising model corresponds to a plaque- 
tte in the original system, and e g means that the corre- 
sponding plaquette operator appears in the decomposition of 

In order to derive Eq. ( 14.1 Oi l, let us define A^p(C|^) 
(A^^(C|g)) to be the number of bonds with ferromagnetically 
(antiferromagnetically) aligned spins in the subset of bonds 
corresponding to g of a given Ising configuration C. De- 
fine as well A^p(e|j) (A^af(C|J)) to be the number of bonds 
with ferromagnetic (antiferromagnetic) spin alignment within 
bonds in the subset complementary to g. Clearly, A^p/af(C?) = 

%AF(e|^)+A^F/AF(e|?). 

We can then rewrite Eq. ( |4.9l l in the Ising language as 

e 




In the following, it is convenient to introduce the conven- 
tion that a bond (ij) belongs to or is inside a partition A of the 
system {(if) E A) if all the spins on the corresponding plaque- 
tte operator belong to A, and the bond does not belong or is 
outside A ({ij) ^ A) otherwise. (Similarly, we will refer to a 
cubic unit cell in or not in A if its six composing plaquettes 
are all in A or not.) 

In conclusion, the numerator of Eq. (14. Il l can be mapped 
onto the partition function of a 3D random-bond Ising model 
on a simple cubic lattice, where the randomness is controlled 
by the choice of g. Again, summation over all possible bound- 
ary conditions is understood. 



Substituting Eq. ( I4.6l l and Eq. ( 14.101 1 into Eq. ( 14. Il l gives 

PiT) = E ^ (4.12) 
geG J ^ > a 

where / = -(1 /2)ln[tanh(PA.B)], Z, = Y^^e^K^A^) is the par- 
tition function of a non-interacting Ising system in a magnetic 
field of reduced strength pA.^, and Zj°'(l) = Zf. 

In the limit of T ^ (p ^ °°), J 0+, all g are equally 
weighed, 

zr(g)=zr(l) VgeG, (4.13) 

and only the states with maximal star magnetization M^{(x) = 
N, i.e., those that are eigenstates of the star operators with 
eigenvalue +1 everywhere, survive: 

gpA,^M,(a) 1 

^ ' ^SM,(a),^- (4.14) 

Such states are of the form g|Oj.), where ^ = 1, . . . , 2^ labels 
the states obtained from |0) by the action of the non-local F 
operators in Eq. ( 12. 3t . Namely, the states jO^;,) are of the form 
r'j'T^^r™^ |0), for all possible choices of mi, 1112, = 0,1. 
The factor 1/23|G| in the above equation appears because 
there are precisely 2^ |G| states with maximal star magneti- 
zation. Thus, one recovers the density matrix of the zero- 
temperature Kitaev model, prepared with equal probability 
across all topological sectors ifisll 

23 

P(7' = 0) = ^tiT^ I ^I0^)(0,|^/. (4.15) 

^ k=i g,g'eG 

In the limit T ^ 00 (p ^ 0), 7 ^ 0°, all g are exponentially 
suppressed except for ^ = 1, while all states a become equally 
weighed. In this case one obtains the mixed-state density ma- 
trix 

^ a 

of a non-interacting Ising model defined on the bonds of a 
simple cubic lattice. 

Clearly from Eq. ( 14.121 ), one expects something to happen 
in the system when the value of the temperature T , i.e., the 
parameter J, is such that the 3D Ising model described by Z]"' 
becomes critical. In order to understand how this relates to the 
presence of topological order at zero temperature, we need to 
proceed with the calculations and compute the von Neumann 
entropy and the topological entropy as a function of tempera- 
ture. 



B. The von Neumann entropy 

Let us consider a generic bipartition of the original system 
§ into subsystems A and S (S = AVJ'E). The von Neumann 
(entanglement) entropy of partition A is given by 

5vN^-Tr[p^lnp^] = -lima„Tr[p:;i], (4.17) 
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where = Tr^p is the reduced density matrix obtained from 
the full density matrix p by tracing out the degrees of freedom 
in subsystem S. Similarly for Sy^, and Sy^ = holds if p 
is a pure state density matrix. 

In order to compute the von Neumann entropy (14.17b from 
the finite-temperature density matrix (14.121 ). we first obtain 
the reduced density matrix of the system using an approach 
similar to the one in Ref. 15, 

Pa^T) = E ^^TTTtE 7 \^A){aA\gA (asl^23l«23) 



pA,,M,(a) 



gee 



l«yi>(ayi|gyi, (4.18) 



where we used the generic tensor decomposition |a) = |ay^) 
l«s)'.? = ^yi®.?S'™dthefactthat {a^\g^\a^) = 1 ifg^ = 
1 3 and zero otherwise. As in the previous section, we denoted 
by Gy^ — {g £ G \ = l,^} the subgroup of G given by all 
operations g that act trivially on S (similarly for G^). 

Notice that a plaquette operator Bp can either act solely on 
spins in partition A (represented in the following by the nota- 
tion p E A), solely on spins in partition "B {p G "B), or simulta- 
neously on spins belonging to A and 23 (which we will refer to 
as boundary plaquette operators, and represent by p G AH). 

Recall from Sec. |lll]that a complete set of generators for 
the subgroup G^^ can be constructed by taking: (i) All pla- 



quette operators that act solely on A, i.e., {Bp \ p S ^1} 

(?v[<J'' = \{Bp I p £ A}\). (ii) All possible (independent) col- 
lective operators constructed from plaquettes in 23 and at the 
boundary, but acting solely on A; as illustrated in Sec. HUl 
the number of such collective operators equals the number rij^ 
of non-contractible loops in subsystem A. And by (iii) ac- 
counting for all constraints given by the independent closed 

(c) 

membranes in A. That is, all 7^}^ cubic unit cells in A, all 
possible (m^ — 1) additional closed membranes if 25 is dis- 
connected, and all independent entire crystal planes inside A 



(m 



(c.p.) 

A 



0, 1,2,3). Again, for all bipartitions of interest in 



our study nij^^'^ 



: and m 



(c.p.) 



3, and for simplicity we will 



restrict to this specific case. 

The cardinalities of the subgroups Gj^ and G^ are thus 
given by 



\Gji\^2^A-^A+"A-('"-B-^) 



■^A 



(4.19a) 



In particular, rtj^ = = 1 in bipartitions 4,5 and zero other- 
wise; and my^ = 2 in bipartition 8, = 2 in bipartition 1, 
and they equal 1 in all other cases. 

Let us then use Eq. ( 14.181 ) to compute the trace of the «-th 
power of Pj^{T): 



MPAiT) 



(«l,yi|^l,yi|a2,yi)(a2 .|g2,.A|a3. 



Al 



K'^n,A\gn.A\'^l,A) 



(4.20) 



Each expectation value above imposes that the two configura- 
tions (i-i^i and a^, / = 1 , . . . , « (with the identification n + 1 = 
1), can be mapped one onto the other over subsystem A, via 
the plaquette flipping operation g^ ^. This is possible only if 
the set , . . . , ^„ e Gyi satisfies the condition HzLi 81 a — ^A' 



i.e., n"=i 81 = ^- Therefore, we can decompose each element 
g, into a product g, = giii^y where |, e G^, / = 1, . . . , n 
with periodic boundary conditions « -I- 1 = 1 (the fact that this 
decomposition is highly non-unique is immaterial to the cal- 
culations below): 



Tr[pM]= L (nSlT 

gy.....g„(^Gj^ \l=l J a^ a„ \/=i 



2-1 111 7 



n«/|0) («L.Alli,yil2,yil«2,yi>(«2,yill2,.4l3,.Al«3,.A)---(««,yill,,,.Ali,yil«i,yi> 
2T{8i) 



1=1 



(n „|3?L^M, (a, ) \ n 
/=i / /=i 



(4.21) 



where we used the fact that the magnetization M^{a.) of state |a) is the same as M^{g(x) of state g\o), for any g e G, to do 
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away with the via relabeling of the states ja^) ^|^,|a,). tion 5(ay^,Pyj^ 

We can further simplify the notation by defining the func- rewritten as 



^yilPyi)' ^^'^ above equation can be 



J 



g,,...,g„GG^ V'=l 

= Z(^)(n)xZ('^)(n). 



\Us^\o)-L n 



/=i 



r 



,a„ \l=i 



1=1 



(4.22) 



Notice that the product n"=i ^ > OC/ 



-l.A/ 



implies 



5(a^^,a„y^), which is therefore redundant and has been 
omitted in the previous equation. In the notation of Eq. (14.221 1. 
it becomes evident that the star (S) contribution, i.e., involving 
only the star coupling constant A-^, and the plaquette iP) con- 
tribution, i.e., involving only the plaquette coupling constant 
A,g, decouple and factorize into two separate terms, Z'^'^'(n) 
and Z(^)(n). In particular, Z(^)(n = 1) = Z'^^\n = 1)^1. 

Using the replica trick, we can compute the von Neumann 
entropy: 



lima„Tr[p:;i] 



lim8„ 



^ -iima„z(^'(«)-iima„z(^)(«) 

= 4'2(^;7'Ab)+4n(-^;7'Aa)- (4.23) 

Thus, from the factorizability in Eq. (14.22b above, it follows 
that the von Neumann entropy has two additive contributions 
from the star and plaquette terms that can then be computed 
separately. 1 10] 

One can check that Eqs. ( 14.22b and (14.23b satisfy indeed the 
r ^ limit discussed in Sec.|lII] as well as the known T ^ °o 
limit (see AppendixlBb. 

Notice that, although in this paper we are concerned with 
3D systems, the derivation is independent of the dimensional- 
ity, and this result holds true for Z2 models in any number of 
dimensions. 

Because the von Neumann entropy is separable as the sum 
of the two independent contributions from star and plaquette 
terms, so is the topological entropy, which is a linear combi- 
nation of the entanglement entropies for the partitions shown 
in Fig. m 



(4.24) 



We now turn to the separate analysis of the two contribu- 
tions. 



(s) 

1. The star contribution Sigpg{T /X^) 

The computation of this contribution is very similar to the 
one in Ref. |4] for the 2D Kitaev model, where the limit Xb 



00 was explicitly considered. 

In order to illustrate this analogy, let us define the following 
entropy differentials: 

= A5(f^(yi;rA^)+A5(S(yi;rAs), 

(4.25) 



and 



A5,opo(r) = 5,„p„(r)-5,„p„(0) 

^topo(^/'^A) + ^'-'topo 



AS[^„,{T/X^) + A^L^pUr/X^), (4.26) 



where 



and 



AS^^'^{A;T/Xs) EE 



S['1{A;T/X^)-S'^I,{A;0) 



AS) 
■*topo 



(0), 



s[^^{A;T/Xj,)-s[^^{A-0) 



Notice that for Xb °°, AS\j^{A;T /Xg) = and 
(p) 

ASi'{T/Xg) = 0. Thus, one obtains that 



AS^yliA;T/XA) = ASy^iA;T] 



(4.27) 
(4.28) 



Moreover, in the limit Xb ^ °° and choosing to work in the 
O' basis, one can show that both the group structure of G and 
the collective operations in Gy[ are very much the same in 2D 
and in 3D. For example, the group G is generated by all but 
one star operators, and the subgroup Gj^ is generated by all 
star operators in A with the addition of all but one collective 
operations that obtain as products of star operators belonging 
to each component of 23 times the ones along the correspond- 
ing boundary. As a result, the topologically non-trivial bipar- 
titions 1 and 4 in 2D correspond to bipartitions 1 and 8 in 3D. 
All calculations generalize straightforwardly to 3D, and one 

(s) (s) 
can derive the expressions for AS'vn and for AS'lopo in a finite 
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system at finite temperature. The actual values for ^(f^ and 

IS) 

■^topo ^re then fixed by matching, say, the known T ^0 limits. 

From the 2D results in Ref. [4], we infer that the star con- 
tribution to the 3D topological entropy is fragile, in the sense 
that it vanishes in the thermodynamic limit at any finite tem- 
perature. Namely, the behavior is singular in that the limits 
of r — > and infinite size do not commute. If the thermody- 
namic limit is taken first. 



T = 
-ln2 r>0. 



(4.29) 



Thus, in the thermodynamic limit, the star contribution to the 
topological entropy evaporates at any infinitesimal tempera- 
ture. 

(The finite temperature and finite size expressions for the 
star contributions to the von Neumann and topological en- 
tropies are shown in Appendix |C]) 



2. The plaquette contribution Sy^{A;T /X^) 

Similarly to the above, one obtains for the plaquette contri- 
butions 



(4.30) 
(4.31) 



Because of the very different nature of the 2D and 3D group 
structures when using the O'^ basis, the computation of the 
plaquette contribution in 3D is not a trivial extension of that 
in 2D, and it thus requires some work. The calculations are 
shown in detail in AppendixiDl while only the results are sum- 
marized here for conciseness and clarity. 

The behavior of A5'topo(r /X^) as a function of temperature, 
in the thermodynamic limit, is 



I T <Tc 
|-ln2 T>Tc, 



(4.32) 



where the critical temperature is associated with a 3D Ising 
transition and can be located at Tc = 1.313346(3)A-b. 



V. DISCUSSION 

We can now put all the pieces together, and argue for the 
persistence of topological order at finite temperatures in the 
3D Kitaev model. Adding the contributions from stars and 
plaquettes, which we have shown to be exactly separable, the 
topological entropy of the system is 



n3D 

^topo 



21n2 T =Q 
(r) = <! In2 < r < Tc 
T > r,.. 



(5.1) 



This is to be contrasted to the 2D case, 



^topo 



(r) 




(5.2) 



where the topological order is fragile, subsiding for any finite 
T (when the thermodynamic limit is taken first). 

In 3D the order survives up to a transition temperature that 
is determined by the coupling constant Xb associated with the 
plaquette degrees of freedom alone. The topological order in 
the system, as measured by the topological entropy, is thus the 
same as in the case where A-a = 0, that is, in a purely classi- 
cal model. In this sense, the order at finite T is classical in 
origin. [22] 

Our results show that the extension of the notion of topo- 
logical order to classical systems applies beyond the hard con- 
strained limit already discussed in Ref. in two dimensions. 
In the 3D example discussed here, the order persists for non- 
infinite couplings A-A, A-B. 

Having obtained the result that topological order in the 3D 
toric code survives thermal fluctuations, in a classical sense, 
up to a finite critical temperature, we now turn to a discussion 
of what this type of order implies. 

At zero temperature, topological sectors can be discerned 
according to the eigenvalues = ± 1 of the loop operators r„, 
where a = 1,2,3, as in Eq. (12.31 1. The eight ground states |/) 
in the different topological sectors can be labeled by integers 
/ = 0, . . . ,2^* - 1 (made up of three bits, / = hhh, la = 0, 1). 

Suppose to prepare, at an initial time t = f,-, a superposition 
of states 



7=0 



(5.3) 



then raise the temperature to some value < T < Tc, and bring 
it back to r = at some time tf. The final T = state will 
again be, assuming thermodynamic equilibrium is reached, 
a superposition of the eight topologically degenerate ground 
states. 

Following the discussion in Section |iy] for temperatures 
below Tc, one can take a winding loop and deform it past ther- 
mal defects, and read off the same eigenvalue of the topolog- 
ical operator as the path is deformed. The information stored 
in all winding loops that do not cross a thermal defect does 
not disappear, as long as there is a way to pass a winding loop 
that avoids defects. Therefore, as long as the system tempera- 
ture is not raised above Tc, upon returning to T = at tf, the 
system should return to the same topological sector that it was 
originally prepared in at time f,-. 

Thus, the state at tf is a superposition 



23-1 

mtf)) = L VFie"" \i) , 



(5.4) 



7=0 



where phases (p/ are accumulated during the thermal cycle. 
These phases, unless locked together by some specific mech- 
anism, shall be randomized by the thermal bath. However, the 
magnitude of the amplitudes remains VpJ, for / = 0, . . . , 2^ — 
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1 , as there have been no transitions between different topolog- 
ical sectors, if the system was never heated above T^. 

Hence, the only (accessible) information preserved under 
the time evolution from f, to f/ is that the relative probabilities 
of find the state in sector / equals pj. The state in Eq. ( I5.4l i re- 
alizes a pbit, or probabilistic bit. lfl3ll It is not a qubit because 
of the thermal dephasing between the states |/). Although still 
a quantum superposition of a sort, in that it has probability pi 
of being in sector /, it cannot be told apart by any type of mea- 
surement from a classical probabilistic system with the same 
probabilities pj. The stability of the system against local mea- 
surements only tells us that the state is not projected onto a 
sector until a non-local measurement is carried out. This effect 
is a non-measurable difference between the state in Eq. ( 15.3b 
and a classical probabilistic state: whether the projection oc- 
curs before (as in the classical state) or after (as in the pbit) 
the measurement is not detectable. 



VI. CONCLUSIONS 

In this paper we have shown that topological order exists 
in the 3D toric code at finite temperatures, up to a critical 
temperature 7^ = 1.313346(3)A-b which is set by one of the 
couplings (that associated to the plaquette terms in the Hamil- 
tonian). This is in sharp contrast to what happens in the 2D 
toric code, where in the thermodynamic limit the order sub- 
sides for any infinitesimal temperature. 

We first presented simple heuristic arguments for this result. 
These arguments are based on the observation that eigenvalues 
of operators defined as products of spin operators along wind- 
ing loops can be used to determine the order even in the pres- 
ence of (thermally activated) local defects, because loops can 
be deformed around such obstacles in 3D, leaving unchanged 
the eigenvalues of such loop operators. This is to be contrasted 
to the 2D case, where one cannot move a loop around a point, 
and thus the eigenvalues of non-local loop operators are un- 
equal on opposite sides of the point defect. 

We subsequently substantiated the heuristic arguments by 
means of an exact calculation of the von Neumann and topo- 
logical entropies in the system as a function of temperature. In 
carrying out this exact calculation, we derived a generic result 
that applies to toric codes defined in any number of spatial di- 
mensions: that the von Neumann entropy is separable as a sum 
of two terms, one associated with stars alone (and a function 
of the dimensionless ratio T /Xa) and another associated with 
plaquettes alone (and a function of the dimensionless ratio 
T /Xb)- The same separability follows naturally for the topo- 
logical entropy, S^^p^iT) = s[^l^{T /X^) + s[^l^{T /Xj^). We 
then showed that, in the thermodynamic limit, the star con- 
(s) 

tribution 5topo(r/A,^) vanishes for any T ^ 0, while the pla- 

(P) 

quette contribution Siq^q{T /Xg) remains constant for T /Xb < 
1.313346(3), and vanishes for temperatures above this scale. 

Because the critical temperature is set by Xb and not Xa, one 
can argue that the topological entropy remains non-zero when 
Xa 0. The resulting Hamiltonian is purely classical, and 
thus one can argue that the nature of the finite T topological 



order must be classical as well. 

Finally, we discussed the nature of the information that can 
be stored robustly in the system because of the topological or- 
der at finite T. We argued that the resilient information stored 
in the 3D system realizes a pbit. 

We end with a note on an interesting situation that should 
occur in systems where both Z2 gauge defect types are con- 
fined. In 3D only one of the defect types is confined, the 
topological entropy drops from 2 In 2 at T = to In 2 for 
< T < Tc, and only the probabilities of being in a given 
topological sector are preserved (magnitude square of the am- 
plitudes, but not the relative phases). If instead both defect 
types are confined, the notion of sectors in both the a" and 
O' basis is retained, and this implies (as discussed briefly in 
Appendix |A]i that, if the system is prepared in a given super- 
position at zero temperature and its temperature is raised and 
again lowered to zero without ever going above Tc, the system 
returns to the same original quantum state (a "boomerang" ef- 
fect). 
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APPENDIX A: THE CONFINED-CONFINED CASE 

In this appendix, we briefly discuss how the nature of 
the topological protection at finite temperature changes when 
both types of thermal defects in a Z2 gauge theory are con- 
fined at low temperature (T < Tc). 

For concreteness and simplicity, let us consider a modifica- 
tion of the 2D toric code, where some ad hoc energy terms 
have been introduced that confine both electric and mag- 
netic thermal defects (without inquiring on the nature of these 
terms. As mentioned in Sec. lIVI this scenario should be real- 
ized in the 4D case without need of any additional term). 

The r = ground state (GS) wavefunction in a given topo- 
logical sector is uniquely specified by the (±) eigenvalues of 
two independent Wilson toric cycles, i.e., winding loop oper- 
ators. In the a' basis, it is sufficient to consider the product 
of all 0^ operators along a horizontal (Tp and a vertical (TJ) 
winding loop, respectively. Similarly, in the basis, using 
loop operators in the dual lattice, and TJ. These loop oper- 
ators satisfy the algebra {71,'J;.} = and {tj,t^} = 0. 

Let us choose to work in the & basis, and define \a,b), 
fl = ±, to be the normalized GS wavefunctions that are also 
eigenvectors of and Tf;, 

\a,b) — a \a,b) 
7l\a,b) =b\a,b). 
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Let us prepare the system in a given superposition of such 
basis states, 

\%n)= L "Va.blaM, (Al) 

a,b=± 

where La,fc=± IVa./;P = 1' ™d consider coupling the system 
to a thermal bath so that the temperature can be varied from 
Tin = 0, via < r < Tc, back to Tfi = 0, as discussed in Sec.llVI 
Trivially, the final state of the system must again be a 
ground state, and therefore it can be written as 

a.b=± 

Moreover, so long as the temperature was never raise beyond 
the deconfining transition at T^, the coupling to the thermal 
bath cannot have transferred any amplitude between any of the 
topological sectors. Hence the following topological quanti- 
ties must be conserved: 

(^,nitf/>i„) = (nitf;.m). (A3) 

For simplicity, consider the case where 

= cos(9/2) 
\|/_.+ = sin(e/2)e''t', 

where 9 G (0, k) and (|) e (— and all others vanish. After a 
little algebra, one can show that the conditions in Eq. (IA3b re- 
quire that the only non-vanishing terms in the final GS wave- 
function are 

\ff+.+ = cos(9/2) 

= sin(e/2)e'*, 



and they satisfy the relations 

cos(e) = cos(9) (A4) 
sin(9)cos((|)) = sin(9)cos(^). (A5) 

That is, 9 = 9 and (|) ±f 

The ambiguity in the sign of (|) is immediately resolved if 
we further require, as expected below Tc, that also the expec- 
tation values of the products / TJ^ and / are conserved, 
leading to the relation 

sin(9)sin((|)) = sin(9) sin(^). (A6) 

Therefore, the quantum topological order in this system is 
fully protected from thermal fluctuations, so long as T < Tc, 
in the sense that the system is bound to come back to the same 
exact initial state upon cooling back to zero temperature. 



APPENDIX B: CHECK AGAINST KNOWN LIMITS 

As a check of the steps leading to Eq. (I4.22l i and ( 14.231 1, let 
us verify that the known limits are indeed recovered. 

For r = (i.e., for 7 = 0) we have that eP^A«.<(«/) /Z, 
KAa,).N/2'\Gl While Zf^g) = Zf\l), yg. In the notation 
introduced below Eq. (I4.14l i. this restricts the summation over 
Ui to states of the form |oC;) = g'l^k)' with g' € G and k = 
1, . . . , 2^ labeling the states obtained from |0) by the action 
of the non local F operators in Eq. (I2.3l l. Namely, the states 
|0j,) are of the form F^'F^^F™ |0), for all possible choices of 
mi,m2,W3 = 0, 1. Eq. (I4.22l i reduces then to 



1 



23"|G|" 



E 

1 

L 



n ^M,ia.),N 
a„ \l=l 



n-l 
1=1 



1*^1 g\,...,g'„eG k,....,k„ 1=1 



1^1 fi'. s;.GG 1=1 



dX' X T^\G\dl-' 



d"A^ X 



n-l 



djj^d^ 



n-l 



r 



(Bl) 



m,)AM+i^k,jA\ ^KmAM+i^)A] 



This in 



where we used the fact that, for the cases of interest, 5 
subsystem A is finite and the non local operators F can 

over , . . . e G can be replaced by an unconstrained sum 



, ^ , ^ turn implies that I g Gm, and the constrained summation 

always be chosen so as to traverse only subsystem i3, / / ^ > 
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mation over g[ & G, g'{, . . . ,g'l, & (where gj^i = g'ig'i+\^ 
for I =2, . . . ,n). Eq. dBlb is indeed the same as in the 2D case 
at zero-temperature. [4] 

In this limit, the von Neumann entropy is given by 



In 



(B2) 



and the topological entropy by ^t^po = 21n2 as discussed in 
Sec.|III](for the full bipartition scheme 1-8). 

For r ^ oo (i.e., for J oo), we have Z f\g) /Zf\\) 
5{g — 1), all a are equally weighed, and Eq. ( |4.22| | reduces to 



n-l 



a^....,a„ l=i 

1 



1 X 



1 X 



•2'iNn 



H-l 



(B3) 



which is indeed the classical entropy of a collection of E^^ free 
Ising spins. The topological entropy vanishes in this limit, 
since the contributions from the different bipartitions cancel 
out exactly (recall that the total number of spins in A for bi- 
partitions 2 and 3 is the same as for bipartitions 1 and 4, and 
similarly for 6,7 and 5,8). 

Notice that, in our chosen factorization scheme in 
Eq. (14.221 1. the plaquette term does not yield any contribu- 
tion to the von Neumann entropy at infinite temperature, while 
at zero temperature the plaquette term contribution equals 
— Int/^, and the star term contribution is —\n{d,g/\G\). 



APPENDIX C: THE STAR CONTRIBUTION 



where E^^ (^3) is the number of a spin degrees of freedom 
in A (23), and E^^ + E^ = 3A^. Here we used the fact that 

5 {cLj involves only subsystem A, hence E^^ spins 

are summed over only once, while there are n independent 
copies of the remaining E^ spins. 
This result leads to 



= -lim3„Tr[p:;i] 



= ln(2' 



= E. In2, 



(B4) 



Here we present the expressions for the star contribution to 
the entropies for finite temperatures and finite system sizes. 
As we argued in the Sec. II V BTI the star contribution to the 
entropies can computed using Eqs. ( 14. 2714. 28b . which relate 
them to entropies evaluated for a hard constrained system 
where — > 00. The calculation in this limit is done most 
conveniently in the O' basis, very much along the lines of the 
calculation carried out for 2D systems in Ref. 4. Paralleling 
the steps of the computation for 2D systems, one obtains for 
the 3D case that 



A5vN(.A;r) 



In cosh i —N 



coshf^fA^-l)) siiihf^(A? 



cosh ( 



coshf^A^) 



-£(i,lnx,) 



cosh(^(A-?^J')) 
cosh (^^a) 



'Y^iyMSi 

i 



sinh(^^(A 
cosh a) 



(CI) 



where Kj^ = -ln[tanh(A.^/r)], JvT^' = +7vri^iB. is the to- 
tal number of star operators acting on the /th component of 
subsystem 25 (either entirely in 25^-, or at its boundary ATi^), 
and 



Notice that only the last two terms in Eq. ( IClb yield a topo- 
logical contribution in our bipartition scheme, since 3sfj'^ — 

•^271 ~ Mi + M/L — Ukewise for bipartitions 5-8. 
Therefore, 



X = cosh f ^ j y = sinh f j (C2a) 
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-(i(2)ini(2) 



1(4 

1=1 

( 

-(x(3)lnx(3) 
+ (xWlniW 

-(xf^'ln^) 
-(x(^)lnx(^) 
+ (x(*^)lnxW 



cosh (^n) 



cosh(^(Af 



cosh 



i=l 

(y' 



(=1 cosh 



cosh(%(A^-?^: 



3A' 



cosh (^A?) 

cosh(^(A^-?^Ji) 



cosh ( 



cosh [^{N-WtL ) 



( 

cosh 




) 


cosh (^^N 
cosh(^(A^-3V 




cosh (^^N 
cosh(§-(A?-]V 




cosh (^^V 





^ coshf^A^j 
, , sinhf^(A?-N^i)) 

^ ' coshf^A^j 

, sinhf^(A?-N<-4)) 
^ ^ coshf^A^j 

, , sinhf^(A?->l'l)) 

^ ^ cosh(^A^) 



sinh(^^(A^ 



^(A^-N^)' 



, Sinn ^uv — JN 
' coshf^A? 



3^' 



:^)iny 



sinh(^(A^-J^;!^)) 
sinh(%(A^-NjJ) 



^ ' cosh ( ^A^ 



(C3) 



where we used the fact that subsystem ® has always one com- 
ponent except for bipartition 1, where it has two components. 

With the expression above for A5topo(r/A-^,A^), one can 
determine the topological entropy contribution from the star 
operators as a function of temperature and system sizes. In 
particular, let us look at two particular limits: that of the zero 
temperature limit taken first, and that of the thermodynamic 
limit taken first. 

For r ^ first, Ka — > 0, and one can easily check that all 
terms in Eq. ( |C3l l vanish, which is expected as the difference 

A5'topo(r/A,^,A^) is, by definition, zero at T = 0. 

Now, when the thermodynamic limit is taken first, i.e.. 



when the sizes A^ and all of 3\f|i);- (for / = 1,2) and N^^^^r ' 
p = 2, . . . , 8 are taken to infinity at fixed Ka, each term in the 
expression in Eq. iC3i gives =Fln2 (with the sign determined 
by whether the partition is added or subtracted). Bipartition 1 
gives —2 In 2 (its contribution is doubled because IS has two 
disconnected components) and it is added to bipartitions 4,5, 
and 8, which give — In 2 each; bipartitions 2,3,6, and 7 are sub- 
tracted, and each of them gives +ln2. Altogether, we obtain 
Is) 

ASiq^q{T /X^,N oo) — — ln2, for any temperature T. There- 
fore, we obtain in the thermodynamic limit the result used in 
Eq. 



One can finally add the zero temperature contributions, to 
obtain 

s['liT/XA) = A5|f)(rA^)-^, (C4) 



and 



s[t{T/K) = A5Lt(rA,)+ln^l444^ 

"23 "3:b "6:b "7S 



-'topo 



topoV 



A5.^o^(rA^)+ln2. 



(C5) 



APPENDIX D: THE PLAQUETTE CONTRIBUTION 



As anticipated in Sec. IIV B 21 the plaquette contribution in 
3D is very different from the 2D case, and we need to carry 
out the calculations explicitly. 

Consider the expression forZ(^), 

- E fri|^l(Oiri^/|0), (Dl) 
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where 



T-tot 



(^)=Lexp 7^ll,7fe)5<5> 
{s,} V (a) J 



is the partition function of the 3D random-bond Ising model 
(summed over all possible boundary conditions), whose ran- 
domness is controlled by g according to Eq. (14. lit . Namely, 
r|,; (g) = ±1 depending on whether the plaquette perpendicu- 
lar to the bond (/j) is flipped in configuration g (r|,y = — 1) or 
not (ri,7 = +1). 

Recall that the group G, and therefore its subgroup G^, is 
defined modulo the identities Dciosed membrane = 1. In the 
language of the randomness realizations {r|,j}, this amounts 
to summing over gauge /nequivalent configurations. In fact, 
any ri,y and rj^ that differ by the product of plaquettes around 
closed surfaces are related by 



3 {5,}. 



(D3) 



Specifically, {S', } corresponds to either of the two spin con- 
figurations that exhibit the closed surfaces in question as their 
only antiferromagnetic boundary (the two configurations are 
related by an overall symmetry). Recall that the product 
of plaquettes belonging to an infinite crystal plane is also an 
allowed gauge transformation, and all possible boundary con- 
ditions (periodic or antiperiodic in each direction) should be 
taken into account when enumerating all configurations {5,}. 
In conclusion, every admits 2^+^ equivalent random- 

ness realizations r|'-^ — y\jjSiSj, labeled by all possible Ising 

configurations {S',}^! (where {S'/Ijlj and {— yield the 
exact same r['-j). 

In the case of a summation over the whole group G, one has 
then the identity 

£ £exp(7£il,-,(g)5,-5, I = 

g^G{S.} \ {ij) ) 

^ E E ( -^Eil'./^'^; I ■ (D4) 

For the subgroup Gj^, the situation is more convoluted. 
First of all, the operators g e G^ correspond to randomness 
realizations {Tl,y(g)} where all the bonds outside A can be 
gauged to assume the value + 1 . Rather than considering all 
the equivalent configurations as for the whole group G, it is 
more convenient to introduce a restricted set of randomness 
realizations {ti|^'} where r||^' is constrained to assume the 
value +1 whenever (/j) ^ A. Notice that we do not constrain 
the bonds inside A, and we are therefore over-counting all the 
gauge equivalent configurations with respect to these bonds. 
The number of equivalent realizations in the restricted sub- 
group can be counted as seen in Sec. |III1 and repeated here- 
after for convenience. All cubic unit cells entirely contained in 
A are independent generators of gauge transformations. Also, 
if A contains crystal planes, there are up to three additional 
generators. Finally, we have one extra generator per con- 
nected component of 23 (i.e., entirely surrounded by JK), but 



for one of them. Thus, the total number of gauge equivalent 
configurations is now 2 



(D2) 

the number of cubic unit cells entirely contained in A, ni 



('"s where again TsT^'' is 

(c.p.) 



is the number of independent crystal planes in A (m 



(c.p.) 



A 



0, 



m^ '''^ = 3 for all cases of interest), and is the number of 
connected components of S. 
As a result, one obtains: 

g'eGJ^{S,} \ {ij) j 



X exp J £ SiSj 



(D5) 



Having done so, the summation over {'n|f''} is now un 



constrained, namely the bond variables r| 



[A) 



±1 are gen- 



erated by freely flipping any of the plaquettes in A, starting 

(A\ 

from the configuration with all rjiy = +1 (which we refer to 

in the following as rj'' = {ilf^}, the ferromagnetic configura- 
tion). Notice that this accounts only for the bipartitions where 
the plaquette operators in A are sufficient to generate the 
whole group Gy^ (bipartitions 1,2,3 and 6,7,8). As discussed 
in SecHni this is not always the case and additional collective 
operations may be needed to generate G^^ (bipartitions 4,5). 
The summation encompasses then all configurations obtained 
by flipping plaquettes in A starting from {ilfy}, and starting 
from the configurations derived from the ferromagnetic one 
via the action of each of the independent collective opera- 
tions. For concreteness, in bipartitions 4 and 5 there is only 
one collective operation in A, illustrated in the bottom panel 
of Fig.|5] In this case, the configurations {ti[j^''} are obtained 
by flipping plaquettes in A starting from the ferromagnetic 
configuration rj", and starting from the configuration with all 
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[A) 



-1, except for those inside the blue thick line in the 



bottom panel of Fig. |5] (i.e., plaquettes in 23 or at the bound- 



ary), where r|) 



[A) 



- 1 . (We will refer to this configuration in 



the following as rj' = {il,^}). If we label {fj'-^' = {fjlj^'}} the 
set of all configurations obtained from the ferromagnetic one 
via the action of the plaquette operators in A alone, the sum- 
mation in Eq. dPSb runs over ri^{fj'-^'} Uri'{fj('^)}, where 
the product of two configurations represents the new config- 
uration with variables given by the site-by-site product of the 

two original variables T|? fi|^' (= fjlj^''), and Tl'^iil^'- 

We can then apply the identity in Eq. ( ID5I ) to simplify our 
expression in Eq. ( IDlb . The condition that a term is non- 
vanishing, namely (0^|g[ ji ■ ■ ■ g'nA\^A) ~ 1' translates into 
the condition that 



(D6) 
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i.e., the product of all 'r\jj (gi), I = 1, n is gauge equiv- 
alent to Tj'' (equivalently ^ ■ ■ ■ s',, a ~ ^^'^y same 
nature of a collective operation in A requires that such opera- 
tion cannot be completed to an identity (a closed membrane) 
by means of plaquette operators in A alone. Therefore the 
above equation holds independently for the collective opera- 
tions and for the fj''^^ configurations. Namely, it imposes that 

the number of collective operations appearing in {f\[f '^}'i^i 
is even, and that 

fl^'h-'\8i)=S.Sj y{ij),3{Si}. (D7) 

1=1 

Trivially, Eq. (ID6b and Eq. (ID7b become equivalent if no col- 



lective operations are present in A. 

Notice that S/Sj = 1 for all (ij) ^ A: all possible {§,} con- 
figurations must be ferromagnetically ordered outside A. If 
ms is the number of connected components in CB, then the 
ferromagnetic order holds across each component separately, 
and from one component to the next the overall sign of the S 
spins may change. An overall sign change of the spins § is 
immaterial, as one can see from Eq. ( ID7I ). and therefore one 
needs to introduce a corresponding factor of 1/2 when sum- 
ming over {Si}. 

Eq. dPlb becomes then 



Z}°'(l)2^V+»'s 



I n 



/=i 



(5,'"} V('7> 



1=1 



.1 



L n 



.n;Lin!/'"=s,s; 



^Z}°'(l)2^V+'"B- 



mi 



(even) 

X E n 

{r|(') };'^ J ('■./> 



1=1 



n 



1=1 



(D8) 



where Lr-m|„ runs over all ntuples {r|''' G ('n''i^')}"=i with an even number of r|' terms. Notice that the summation 



V '=1 J 



(D9) 



where Z„ (^{yfflj'^f^^y'}; can be interpreted as the partition function of an Ising chain of degrees of freedom 
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{^\f''^}'i=i in a random field of local strength Jylfj'sf^sf, 
and subject to the condition that the product of all Ising spins 

Yl"=i'Hij^''^ equals SiSj. By means of the change of variables 
^(yi,/) _ ^{Aj)^{Aj+i)^ jj^j^ becomes the partition function 
of a nearest-neighbor Ising chain with periodic or antiperiodic 
boundary conditions (BC) depending on the sign of S^Sj — ±1 

(i.e., m\- 



2j„ — 



z„{{ji?sf^sf};s^s;) 



BC=SiSj 



This in turn can be computed exactly. 



2Z„ = (2coshy)"- 
= (2coshy)" 



StSj 



(DIO) 



(2sinh7)" 



We also used the fact that ff|'.'' = +1 if (ij) e yi by construc- 
tion. (Notice that this convenient choice does not introduce 
any limitations. In general, the number of times when a —1 
appears in the / = 1 , . . . , n sequence of rf|j' values must be 

even, and therefore Il/Li = +1, V/,7). 



11)^(1) 



SiSj 



(2sinh7)" . 



(DID 




For convenience of notation, let us consider the following 
change of summation variables 



(D12) 



so that we can write Zn — \ e'^" e^"^'^i, with A„ and B„ defined 
as 

gA„+B„ ^ (2coshy)" + (2sinh7)" (D13a) 
gA„-B„ ^ (2coshy)" - (2sinh7)" (D13b) 



Given that HzLi'^f'' = il. for all sites / whose adjacent 
bonds (//') are solely in A, the summation over {Sj = ±1} and 
the summation over {9, = ±1} are unconstrained. The case 
is different for the sites / that have an adjacent bond not in 
A. The correlation across such bond is in fact ferromagnetic 
by construction, and, if 25 has only one connected component, 
the spin 5, has the same sign as all other spins not entirely 
surrounded by bonds in A. Consequently, all the boundary 
spins Si have the same sign, and the values of the associated 

spins 9; are determined uniquely by the product 11"= 1 s!f\ If 
is the number of connected components in !B, then the 
ferromagnetic order holds across each component separately, 
and from one component to the next the overall sign of the 
S spins may change. This is accounted for by summing over 
boundary sign variables ^'^ = ±1, j^ws, assigned to 

each boundary 3^ defined as the set of sites that have adjacent 
bonds both in A and in the rth component of A. 
In the end, Eq. (ID8b becomes 



1 



^Z)°'(l)2^V+'"s-i 
1 



(even) 



'] I .1 



zr(i)2' 



1 1 (even) 

^ E L n L n 



{{sf}]" {e.}{';>e-^ {r|(')}'/^i('»^^ 



expl/Ln^P^? 



1=1 



kr=±l}:i? '-=1 '^3. V '=1 J 



1 



^■"a ^" 1 



(even) 



Zf (1)2^ 



2^i^' 2 



^ E E n ^"""'''^ E n 

{{sf}}" {e'}<'^'>e-^ {tt(')}'/^,('7>0^ 



exp UEn^P^f 

V 1=1 , 



X E uu^^.nsf'-^r 



(D14) 



Notice that ri|'' = +1 if the plaquette (/;) does not belong to the collective operation, and that whenever (ij) belongs to 
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the collective operation the value of r||j' = ±1 is the same for 
all (ij) . (We restrict here for simplicity to the case where there 
is at most one collective operation in A. In order to extend to 
the general case one needs to repeat the derivation for each 
collective operation separately.) 

Notice also that the sum over sf^ that are entirely sur- 
rounded by bonds in A is unconstrained, and it contributes a 

trivial factor 2^ " to the sum over the remaining spins. In the 
following, we use this simplification and all summations over 
5,-'' are intended as constrained only to the remaining spins 
(for convenience, we do not increase the already complex no- 
tation). 

Let us focus on the boundary condition 



g^=±l}-B r=l,Ga,- V 1=1 



{') 



(D15) 



where [ J stands for the integer part of its argument. In other 
words, the sum 9/ + YJi=i must equal n + q (mod4), or 



e, + £ Sf^ - (n + q,-) = (mod 4) . 
1=1 



Using the function 



^ k=0 ^ ^ ' 

_ J 1 if X = (mod4) 
"jo if x= 1,2,3 (mod4) 



(D16) 



(D17) 



Given that the 9 and the S spins can assume only the values 
± 1 , then the quantity 9; + YIi= \ ^ can only assume the values 
n + l,n — l,n — 3,..., — (n — 1),— (« + l). ll23ll In particular, 

the product QiYl'i=\Sf^ is positive whenever said summation 
equals «+l,n — 3,« — 7..., and it is negative otherwise. We 
can therefore rewrite the delta function in the above equation 
as 



V 1=1 



[in+qr)/2\ ( n \ 

p=Q \ 1=1 J 



we can finally write the delta function as 



h\%,Wsf=q, 



1=1 



llexp 



i'^kiU + psf^-in + q,-) 
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Substituting into Eq. (ID14l i. we obtain 



zr(i)2^ 



3v''''a„ , , 1 (even) 



i'^kAQ, + Y,sl'^-{n + qr) 



1=1 



1 



" e^A^^ 1 1 



-1 .T-J, H L 



(even) 

X E n L L 



exp u L nij'^r^^r 



exp / L L^sPsf 

. V {ij)Ml=2 , 



exp< 



/g3 



\Y.Y.^\^^^s'^P-\-qr 
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where 8 and H-^ are, respectively, the full set and the total both in A and outside A. In the language introduced earlier, 
number of boundary sites, i.e., sites that have adjacent bonds _ _|_ jvj-W _ ^ j^(<;) ^ therefore Jsf— = 

A A 'd'' j\ 
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^3 +^a- 



Note that the last line in Eq. ( ID 191 ) does not depend on the 
5'') or 9 spins. If we introduce the partition functions 



{5,} 



(ij)iA 



ied 



{s,} 



we can carry out the summation over the even number of col- 
lective operations {tT^''}"=i explicitly, and arrive at 



= ( ^ , 



A ^" 1 



2 



£ exp Bn £ 9,9^ ) exp 

{e,} V (r7>6yi 

+ 



''=1 /g3,- 



exp 7 E 5lXM-exp y £ 



(4; 
(4; 



,1 c(l)c(l) 



(D20) 



We are finally in the position to take the derivative with respect to n, and to compute the von Neumann entropy of the bipartition 



s'y^^{A;T /Xs) = -lima„Z(^)(«) 



M— >1 



1 



(l)c(l) 



exp U L ^r^^) 



liJ^ ^ / [{9/} V ('i>G^ / 



Um 3„ 



exp 7 E S\'^sf^ 



i E i E exp 



z)°'(i)2'"s-i y 2 



£exp Bi £ 9,9, 



E ^ E i E e^p 



^ {9,} \ (!7>GA 



r=l iear 



exp| / E 5!"5yM+exp y E ^1^^ 

expfy i: ^p'^y^Vexpfy i: ^^^^^A 



(D21a) 



(D21b) 



(D21c) 



r 



The summation over {^,} can be carried out explicitly Eq. ( ID21l i. This leads to a delta function that identifies 
both in the first (|D21a|l and second (|D21b|l contribution to q. ^ ^^5^)^ / e 3^ and r = 1, . . . One can verify that 



22 



the factor qr is actually immaterial, and the 9 and 5^'^ terms summation over {qr = ±1}'^[ becomes then trivial, yielding 
in the above equation can be gathered into a single partition an overall factor 2™* . 



function 



{e,} \ {ii)eA J r.(l), V (ij)iA J 



= £exp[y£5,5j^zr(l), 



(D22) 



where we used the fact that Bi = J (see Eq. ( ID26b below). The 



In the third contribution jD21cl ). each summation over qr = 
±1 yields a factor 2cos(7i5^,ga^J:,/2), which vanishes unless 
Y,ki is even. Thus, we can constrain the summation over 
{ki = 0, . . . ,3},ga,. to satisfy this condition, and we can drop 
the terms exp [/^ Liea, ki{l— qr)\ , since 1 — q'r is even and the 
term is identically one. The summation over {qr — ±l}"j°j 
becomes again trivial. In particular. 



exp 



exp^/-^^, (9,-1) +(5,'''-! 



(D23) 



for the same reasoning, and we can write the 9 and 5*^^^ terms (The labeling !B instead of A is used here as a reminder that 

the summation over {9, } includes both spins surrounded only 
by bonds in A, and spins on the boundary 3. Therefore, the 



in a more compact form using the definition of •Z^^j'^, and in- 
troducing the notation 



{e,} 



J E 9,9,+/- ^^,(9,-1) 

{ii)eA ied 



(D24) 



total number of 9 spins is:Mj' =k^>+:m3.) 



These considerations allow us to simplify Eq. ( ID21b to 



S\;>,{A-T /Is) = -lima,, 



1 



(even) Z 



,Z)°'(l)2"'B-i 



n-1 



2^A 

£exp(B„ SiSj+J ^ S^Sj 

,{S,} V ('7>G-A {iJ)M y 



(p) 



1 



lim 3„ 



1 



(D25a) 



(D25b) 



{ (4a +4J ) 1" (4a +4a ) + (4f -4} ) 1" (4f -4a ) } • ^^^^so 



In order to proceed further, let us first study some of the 
terms in Eq. ( ID25b separately. From Eqs. ( ID 131 ) we have that 



A„ = -ln{ [(2coshy)" + (2sinhy)"] 
X [(2coshy)"-(2sinh7)"]} 
= iln|(2coshy)2"-(2sinh7)^"| (D26a) 



B„ 


1 

= -In 
2 


(2coshy)"4 


(2sinh/) 


(2 cosh/)" - 


(2sinh/) 


Ai 


= ln2 






Bx 


^ lln 
2 


1 + tanhy 
1 - tanhy 


J 



(D26b) 
(D26c) 
(D26d) 
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dn 



d 

-l-Bn 

dn 



= ln2 + cosh^y In (cosh7) — sinh^7 In (sinh7) 



= sinhy coshy In - 



sinh/ 
cosh 7 



(D26e) 
(D26f) 



Notice that 



In 2 for J 



1/2 + {e-^-') for 7 oo, and that ^B„ I ^ 



0, i;A„ 



n=l 



-7 + 
for 7 ^ 0, 



-l/2 + 0{e-^-') for7^oo. 



We can also carry out the derivative in Eq. ( ID25bb : 



J 



lim d„ 



^exp B„ SiSj+J S.Sj 



US,} 



(ii)iA 



d 

-T^n 

dn 



Y\ L 5,5/ exp Bi Y SiSj + J L S^Sj 



= 1 {S,} \(ij)eA 

sinh 7 



sinh 7 cosh 7 In ■ 



cosh 7 



Y L 5,5, exp 7^5,5, 



[Si] \(ij)eA 



{ij) 



■^A)zf'{l) 



sinh7 mt/ X 

sinh7cosh71n — — (fi^) Zy (1) 
cosh7 ^■1 

(l(,7)g^^.^;)exp(7I(,^.^5,-5, 
Zf(l) 



(D27a) 



(D27b) 



where Ej^ is the extensive energy of the bonds in A (in units of 7), in the Ising model described by the equilibrium partition 
function Z}°'(1). 

The last calculation we still need is 



dn 



ln2 + cosh^7 In (cosh7) — sinh^7 In (sinh7) 



Combining all the results in Eqs. ( ID26l l. (D21\ and ( ID28I ). Eq. (D25\ reduces to 



(D28) 



S^^''^{A-T /Xb) = ln( 2"'«-' ) +lnZr(l)-:M; 



sinh 7 

sinh7cosh71n : — (E 



ln2 + cosh^7 In (cosh7) — sinh^7 In (sinh7) 



cosh7 ^'-^'zfd) 



1 1 



Z)°'(l)4^^ 



(D29a) 
(D29b) 



+ _ yA- 
} ^{k} 



(D29c) 



Recall that Y.h is even, and therefore L^;5, is also even, and both Z^^t^ and zt^^ can be rewritten as 
irrespective of the values of the spins {5, = ± 1 } . In particular. 



exp 



n .-'■?^'cos^fc, + /.-'1^' (sin^fc,)5, 



■ [5*:, even + 5,- hk^ odd] , 

ig3 



(D30) 



yA,\- 



jA,+ 



= Z 



A, odd \ / > 

{5,} V iG3 / \ (/;>Gyi / 
//t,odd \ 

'^'^ { n 

A-,odd \ / > 

{S,} V /g3 / V {ij)<^A / 
/ kf odd 



/g3 / 



(D31) 



(D32) 
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where Z^'+ = I{,_.}exp(7I^..j^^5,5,.) and Z^-+^ = 

L{s,}^w{j'L(,j)^ASiSj)- Similarly for Z^"'^;^ and Z^^^;^. 
Thus, all these quantities can be interpreted as correlation 
functions of boundary spins located at the odd entries of the 
set {ki} times a partition function. Note that the constraint 



Y^iedr^i even, Vr, requires that the number of such odd 
entries is also even separately on each boundary component 
r= l,...,ms. 

If we are interested in computing the topological entropy 
of the system, it is convenient to decompose the last term in 
Eq. (Ipj^l l so that 



In 2 



-lnzr(l)-?^if^ 



ln2 + cosh^y In (cosh/) — sinh^7 In (sinh7) 



sinhy 

smh/cosh/ln— ^(£^>^^.,1^ 
4^a \ zr(l) 



J (even) zf.'t Zf' 



zr(i) 2 



In 1 
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3} 

7^1,+ 



In 1 



(D33a) 
(D33b) 

(D33c) 



(D33d) 



The result in Eq. ( ID33l l holds for ny^ = 1 (i.e., there is only 
one collective operation in A). In order to compute the topo- 
logical entropy of the system with the bipartition scheme in 
SecHni we also need to consider the case where = 0. Re- 

I 



and 



Notice that Eq. ( ID35l l differs from Eq. ( ID33l l only in that it 
lacks its last contribution ( |D33dt . 

We can finally compute the plaquette contribution to the 



peating the derivation above, from Eq. ( ID19I ) to Eq. ( ID33I ), in 
the absence of collective operations leads rather straightfor- 
wardly to the result that 



(D34) 

(D35a) 
(D35b) 

(D35c) 



topological entropy S\q^q{T /Xg), using the full bipartition 
scheme. All the terms that do not carry a topological con- 



z(^)(„) = ( L_ 



-1 



n ■|\r''''4 , 

e^A 1 



E i E Eexpffi« E e-e.) 



E 



' r=liedr 



l)c(l) 



exp / E Sm 



j-1 



S^^^{A-T /-kg) 



= ln(2"'«-i ) +lnZ}°'(l)-3^; 



sinhy 

smhy coshy In (E 



(p) 

A 



ln2 + cosh^7 In (cosh7) — sinh^7 In (sinh7) 



coshy ^^-^^zrd) 



1 (even) y'S>-+ yJ^& 
4^3 \ Zf{\) {^'} 
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tribution cancel. Namely, as discussed in Sec. Hill 
and on similar grounds 



Likewise for bipartitions 5-8. Recall also that = 1 and 
nji~Q for all bipartitions, except bipartitions 4 and 5 (which 
(D36a) have = 1 and nyi = 1), and bipartition 1, (which has — 
2 and ha — 0). Using Eq. ( ID33I ) and Eq. ( ID35I ) accordingly, 
we obtain 



lAI 



(D36b) 



J 



In (2"™i'-B+"'2s+« 

1 (even) fzjf^+zlf. 



4^a \ Zr(l) 2 
(partitions 5 ... 8). 



{k,} 



y2'B.+ y2A.+ y3'B.+ y3A.^ 



zr(i) 



74.A,- 
AA,+ 



In 1 + 



y4yi,H 



zr(i) 



InZ 



y^AA 



In 1 



74®, 



^4^1, 



3^1,+ I {kj} {kj} 
{k,] + 2f (1) 



^4^1.- 
,471,4 



InZ' 



471,4 



(D37) 



Using the fact that — ~ «^ 



3S " 



Ms 



1 , that m 



5S 



7S 



'8S 



= 0, and that 'z!]j^,^ (J) = Z^^^"*" (7) since 



bipartitions 4 and 5 are in fact identical, one arrives to the 
result 



CpUrAfi) = -ln2 



-j^ (even) 

4 3 X, 
{k.},=\ 

2 (even) 

4 3 X, 
{k.},=] 



yl'B.+ ylA.+ 

tiM_li^lnz!f;^ 

7tot/t\ [kij 



zr(i) 



753.+ ^571.+ 

75/1,4 

zf(i) ^''^ 



72B.+ -72/1,+ 

7totti\ {ki} 



zr(i) 



76S.+ y6A.+ 



zr(i) 



^ (even) Z^f .+ Zf^-^ 
^ {kit 

4 3 TM-. 



zr(i) 



733,+ 73.A,+ 



zr(i) 



yTS.+ ylA.+ 




74/1,- 



^{kj} ^{kj} ,„^4A,H 



783.+ 78.^,+ 



In 1 - ^ 



^4.4,- 
{k,} 



-AAA 
^{ki\ 



^{k,} 



(D38) 



This expression can be cast in a more useful way by notic- for each of the partitions p = 1, ... ,8. This is because the 
ing the following. Factors like expectation values of the products of spins is always non- 
negative because the interactions are ferromagnetic (this can 

1 Z'f^i^ Zfj^{^ be shown explicitly in a high temperature expansion, for ex- 
CPn, I = — ' J ' — (D39) ample). Recall that the set {kj} contains always an even num- 

4 3" ^/ ^ ^ ^ berofoddyti's. 

2 j^p'h,^ 2^pAA 




4^3P Zf{\) 

> 0, Moreover, one can check that 
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(even) 



^ > {e,}{s,} V (u>GP^ / V {mp-A / 

^expjy ^ e,e^.)exp(y ^ 5,5 J 

{S,} V {ij)epA J \ {ij)^pA ) 



Y (even) 



4^3" 



2 (even) 
4^3" 



/^EA:,(e,+5,-2) 



= ^I:L^^p(^ L e,0Jexp(7 Y. ^--5. U L i E exp /|£^,(e,+5,-l-^) 

V'''J;=1 



{S,} V {ii)eA J \ (ij)<^A J ^c,=±l 

^xzr(i) = i, 



zr( 



A' {e,}f- 



(D40) 



and thus the J'j'j. } > are probabiHty weights. 
Similarly, we can define a probability 

= f ^2 j,3 j,6 j,7\ 

> 0, (D42) 
where the {ki} are defined on the total boundary of the added 



partitions, and we used the fact that partitions 1,4,5,8 and 
2,3,6,7 have exactly the same total boundary. We can then 
define averages with respect to this measure, 

(even) 

L %}(•••), (D43) 

and Eq. (ID38l l reduces to 




+ In m 1^ (D44a) 



(D44b) 



We can finally analyze this expression as a function of 
temperature. Recall that J = — (l/2)ln[tanh(PA-g)], so that 
7-^0 when J — > 0, and the disordered Ising phase occurs 
for T <Tc 1.313346(3)A-B. Below the Ising transition at 
J = Jc — 0.2216544(3), one can use a high-temperature loop 

expansion to estimate the ratio of Z^^j over Z^^j^. 

The high-temperature expansion contains either closed 
loops, or open strings that terminate at the boundary, because 



an Si is inserted for each site / where ki is odd. The corre- 
sponding expansions for Zi^j over differ only by loop 
terms that intersect the twist surface (generated by the collec- 
tive operation in Fig.|5]bottom) an odd number of times. These 
terms appear indeed with opposite sign in the two expansions. 
This can be achieved only by closed loops that wind around 
the donut shape, and by open strings that connect boundary 
spins Si among those identified by the set of odd A:,'s (see 
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Fig.©. 



twist surface 




FIG. 8: (Color online) - Qualitative examples of terms in the loop 
expansion that appear with different signs in Zj^j and Z|^j'^: 
closed loops that wind around the donut shape, and open strings that 
connect boundary spins 5/ (which appear in the high-temperature ex- 
pansion whenever the corresponding k, is odd). 

In the high temperature limit, long loops are exponentially 
suppressed and we can safely neglect the winding loop contri- 
butions when the size of the partition is taken to infinity. Sim- 
ilarly, out of all possible ways of connecting boundary spins 
in the ki odd set, only 'short' strings between spins 'close' to 
the twist surface need be considered, as illustrated in Fig.|9] 




FIG. 9: (Color online) - Schematic, projected illustration of open 
strings between boundary spins. The location of the spins 1 , 2, .... 8 
are given by the sites where k, is odd (recall that their total number 
must be even). One can verify that the parity of the number of inter- 
sections with the twist surface is fixed by the choice of the locations 
1 , 2, . . . , 8, up to exponentially small corrections such as the red dot- 
ted string in the figure, which vanish in the thermodynamic limit of 
A'g oo. For example, consider the change upon reconnecting spin 
5, . . . , 8 via the dashed lines instead of the solid lines. (Notice that 
the case where, say, the points 1 , . . . , 4 are uniformly distributed on 
the boundary is exponentially suppressed by the probability ^Pji,}.) 

For ki points near the twist surface, rearranging the way 
that points are paired does not change the parity of the num- 



ber of crossings of the twist surface. This is illustrated in 
Fig.|9] where reconnecting spins 5, ... ,8 via the dashed lines 
instead of the solid lines give instead of 2 crossings, thus not 
changing the parity. Now, a reconnection that changes the par- 
ity involves drawing long strings. Below the Ising transition, 
the probability y{ki} keeps the points with odd ki confined in 
pairs; thus there are ways to connect them together with short 
strings. But changing the parity of the intersections requires 
re-matching them in such a way that connections with sites far 
away are made, and the total length of these strings is of or- 
der the system size. This is illustrated in Fig.|9l for example, 
reconnecting spins 1, ... ,4 requires strings whose total length 
spans the system size. 

Therefore, one can verify that all the loop terms corre- 
sponding to a given choice of fc,'s have the same parity in the 
number of intersections to the twist surface, up to corrections 
that are exponentially small in the size of the bipartition. As 

a result, the ratio Z'^'^j^ /Z'^'^j^ tends to ±1 in the thermody- 
namic limit of — > oo, and the sign is purely determined by 
the choice of ki. 

Eq. ( ID44bl ) is clearly symmetric under the change 
Z^^j^/Z^^j^ — Z|^j^/Z|^j^, and we finally arrive at the 
result that at low temperature T < Tc, the term in Eq. ( |D44b| ) 
gives 21n2. 

In the Ising ordered phase (T > Tc here), on the other hand, 
the ratio Z^'^'^/Z^'^'+ ^ in the thermodynamic limit, be- 
cause of the energy cost associated with the twist in boundary 
condition (domain wall) in the '— ' partition. Hence, in this 
case the term in Eq. ( |D44b| ) gives 0. 

A similar reasoning gives that the ratios entering Eq. ( |D44a| ) 
are equal to 1 in the thermodynamic limit, and corrections 
appear only as the coiTelation length becomes of the order of 
the size of the bipartitions, i.e., infinite in the thermodynamic 
limit. Thus, in the low temperature phase, Eq. (ID44al ) gives 
Inl =Oforr < r,.. 

On the other hand, for T > Tc, the partitions order ferro- 
magnetically, and one must account for the fact that partition 
lA has two disconnected components, and therefore these two 
components can order in two ways relative to one another, 
giving a factor of 2 in the ratio appearing in Eq. (ID44al l. and 
hence this terms gives a contribution ln2. 

Putting it all together, we obtain that 

and AslZiT/K) = <o(7^/^b) - ^topUo) is given by 
Eq. ( I43a . 
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